      SUBROUTINE RADA27(X,V,TI,TF,XL,LL,NV,NI,NF,NS,NCLASS,NOR)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!  E.EVERHART, PHISICS AND ASTRONOMY DEPARTMENT, UNIVERSITY OF DENVER 74
!  DENVER,COLORAD 80210. PHONE (303)-753-2238 OR 753-2362
!  PROGRAM ADAPTED FOR BESM-6 BY S.TARASEVICH, DEPT. OF ALGORITHMIZATION
!  INSTITUTE OF THEORETICAL ASTRONOMY, LENINGRAD 192187 USSR.  1975
!  INTEGRATOR FOR ORDERS 7,11,15,19,23,27.  NOR IS THE ORDER.
!  DOUBLE PRECISION VERSION FOR BESM-6  OF INTEGRATOR RADAU
!  NV IS THE NUMBER OF DEPENDENT VARIABLES
!  NCLASS IS 1 FOR 1ST-ORDER DIFF EQ,AND 2 FOR 2ND ORDER DIFF EQ.
!  IF FIRST DERIVATIVES ARE NOT PRESENT (CLASS IIS),THEN USE NCLASS=-2.
!  X(NV) IS THE INITIAL POSITION VECTOR.  IT RETURNS AS THE FINAL VALUE
!  V(NV) IS THE INITIAL VELOCITY VECTOP. IT RETURNS AS THE FINAL VALUE
!  IN THE CASE NCLASS IS UNITY,THEN V IS SIMPLY ZEROED
!  TI IS  INITIAL TIME, TF IS FINAL TIME, NF IS NUMBER OF FUNCTION EVALU
!  NS IS THE NUMBER OF SEQUENCES.
!  PROGRAM SET UP FOR ANY NUMBER OF SIMULTANEOUS EQUATIONS
!  LL CONTROLS SEQLENCE SIZE. THUS SS=10.**(-LL) IS DESIRED SIZE OF A TE
!  AS IN CONTROL SYSTEM I.
!  IF LL.LT.0, THEN XL IS THE SPECIFIED CONSTANT SEQUENCE SIZE
!  WILL INTEGRATE IN A NEGATIVE DIRECTION IF TF.LT.TI
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      DOUBLE PRECISION X(NV),V(NV),F1(36),FJ(36),C(78),D(78),R(78),
     a XI(78),Y(36),Z(36),
     a B(13,36),BE(13,36),H(14),W(13),U(13),BT(13,36),HH(48)
      DOUBLE PRECISION ZERO, ONE, TI, TF, TM
      DOUBLE PRECISION T,XL,S,Q,TEMP,RES,T2,W1,TP,SR
      DOUBLE PRECISION VAL,TVAL,HSUM,SM,BDUBL,PW,TDIF,DIR,SS
      COMMON /F1/F1 /FJ/FJ /Y/Y /Z/Z 
      COMMON /BE/BE /BT/BT /B/B 
      DIMENSION MC(12),NW(14),NXI(78)
      LOGICAL J2,NPQ,NSF,NPER,NCL,NES
      DATA NW/0,0,1,3,6,10,15,21,28,36,45,55,66,78/
      DATA MC/1,13,24,34,43,51,58,64,69,73,76,78/
      DATA ZERO,ONE/0.D0,1.D0/
      DATA NXI/2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 3, 6, 10, 15,
     c 21, 28, 36, 45, 55, 66,
     c 78, 4, 10, 20, 35, 56, 84, 120, 165, 220, 286, 5, 15, 35, 70,
     c 126, 210, 330, 495, 715, 6, 21, 56, 126, 252, 462, 792, 1287,
     c 7, 28, 84, 210, 462, 924, 1716, 8, 36,
     c120,330,792,1716,9,45,165,495,1287,10,55,220,715,11,66,286,12,
     c78,13/

      DATA  HH/.212340538239152D0,.590533135559265D0,.911412040487296D0,
     c.0985350857988264D0,.3045357266463639D0,.5620251897526139D0,
     c.8019865821263918D0,.9601901429485313D0,.0562625605369221D0,
     c.1802406917368924D0,.3526247171131696D0,.5471536263305554D0,
     c.7342101772154105D0,.8853209468390958D0,.9775206135612875D0,
     c.0362578128832095D0,.1180789787899987D0,.2371769848149604D0,
     c.3818827653047060D0,.5380295989189891D0,.6903324200723622D0,
     c.8238833438370047D0,.9256126102908040D0,.9855875903511235D0,
     c.0252736203975203D0,.0830416134474051D0,.1691751003771814D0,
     c.2777967151090321D0,.4015027202328608D0,.5318623869104160D0,
     c.6599918420853348D0,.7771593929561621D0,.8753807748555569D0,
     c.9479645488728194D0,.9899817195383196D0,.0186103650109879D0,
     c.0614755408992690D0,.1263051786933106D0,.2098429717265625D0,
     c.3078989982803983D0,.4155560359786595D0,.5274156139958823D0,
     c.6378686027177612D0,.7413764592942375D0,.8327489886084423D0,
     c.9074047753009974D0,.9616018612603216D0,.9926353489739107D0/
      KD=(NOR-3)/2
      KD2=KD/2
      KE=KD+1
      KF=KD+2
      PW=ONE/(KD+3)
      NPER=.FALSE.
      NSF=.FALSE.
      NCL=NCLASS.EQ.1
      NPQ=NCLASS.LT.2
      SR=1.5
      IF(NV.EQ.1) SR=1.2
      NES=LL.LT.0
      TDIF=TF-TI
      DIR=TDIF/DABS(TDIF)
      IF(NES) XL=DABS(XL)*DIR
      NCLASS=IABS(NCLASS)
      LA=KD2*KD2-1
      DO 14 N=2,KF
      LA=LA+1
      H(N)=HH(LA)
      W(N-1)=ONE/(N+N**2*(NCLASS-1))
  14  U(N-1)=N+1
      DO 22 K=1,NV
      IF(NCL) V(K)=ZERO
      DO 22 L=1,KE
      BT(L,K)=ZERO
  22  B(L,K)=ZERO
      W1=ONE/NCLASS
      DO 939 J=1,KD
      M=MC(J)
      JD=J+1
      DO 939 L=JD,KE
      XI(M)=NXI(M)*W(J)/W(L)
 939  M=M+1
      C(1)=-H(2)*W(1)
      D(1)=H(2)/W(2)
      R(1)=ONE/(H(3)-H(2))
      LA=1
      LC=1
      DO 73 K=3,KE
      LB=LA
      LA=LC+1
      LC=NW(K+1)
      JD=LC-LA
      C(LA)=-H(K)*C(LB)
      C(LC)=(C(LA-1)/W(JD)-H(K))*W(JD+1)
      D(LA)=H(2)*D(LB)*W(K-1)/W(K)
      D(LC)=(D(LA-1)*W(K-1)+H(K))/W(K)
      R(LA)=ONE/(H(K+1)-H(2))
      R(LC)=ONE/(H(K+1)-H(K))
      IF(K.EQ.3) GO TO 73
      DO 72 L=4,K
      LD=LA+L-3
      LE=LB+L-4
      JDM=LD-LA
      C(LD)=W(JDM+1)*C(LE)/W(JDM)-H(K)*C(LE+1)
      D(LD)=(D(LE)+H(L-1)*D(LE+1))*W(K-1)/W(K)
  72  R(LD)=ONE/(H(K+1)-H(L-1))
  73  CONTINUE
      SS=(1.D1)**(-LL)
      NL=NI+30
!  SET IN A RESONABLE ESTIMATE TO T BASED ON EXPERIENCE.  (SAME SIGN AS
      TP=((NOR/1.1D1)*(0.5D0)**(0.4D0*LL))*DIR*5.D-1
      IF(NES) TP=XL
      IF(TP/TDIF.GT.0.5D0) TP=0.5D0*TDIF
      NF=0
      NCOUNT=0
4000  NS=0
      TM=TI
      SM=1.D4
      CALL FORCE(X,V,TM,F1)
      NF=NF+1
!  LOOP 58 FINDS THE BETA-VALUES FROM THE CORRECTED B-VALUES, USING D-C
 722  DO 58 K=1,NV
      BE(KE,K)=B(KE,K)/W(KE)
      DO 58 J=1,KD
      JD=J+1
      BE(J,K)=B(J,K)/W(J)
      DO 58 L=JD,KE
      N=NW(L)+J
  58  BE(J,K)=BE(J,K)+D(N)*B(L,K)
      T=TP
      TVAL=DABS(T)
      T2=T**NCLASS
      DO 175 M=1,NL
      J2=.TRUE.
      DO 174 J=2,KF
      JD=J-1
      LA=NW(JD)
      JDM=J-2
      S=H(J)
      Q=S**(NCLASS-1)
      IF(NPQ) GO TO 5100
      DO 1300 K=1,NV
      RES=B(KE,K)
      TEMP=RES*U(KE)
      DO 7340 L=1,KD
      JR=KE-L
      RES=B(JR,K)+S*RES
7340  TEMP=B(JR,K)*U(JR)+S*TEMP
      Y(K)=X(K)+Q*(T*V(K)+T2*S*(F1(K)*W1+S*RES))
1300  Z(K)=V(K)+S*T*(F1(K)+S*TEMP)
      GO TO 5200
5100  DO 1400 K=1,NV
      RES=B(KE,K)
      DO 2340 L=1,KD
      JR=KE-L
2340  RES=B(JR,K)+S*RES
1400  Y(K)=X(K)+Q*(T*V(K)+T2*S*(F1(K)*W1+S*RES))
5200  CONTINUE
      CALL FORCE(Y,Z,TM+S*T,FJ)
      NF=NF+1
      IF(J2) GO TO 702
      DO 471 K=1,NV
      TEMP=BE(JD,K)
      RES=(FJ(K)-F1(K))/S
      N=LA
      DO 134 L=1,JDM
      N=N+1
 134  RES=(RES-BE(L,K))*R(N)
      BE(JD,K)=RES
      TEMP=RES-TEMP
      B(JD,K)=B(JD,K)+TEMP*W(JD)
      N=LA
      DO 471 L=1,JDM
      N=N+1
 471  B(L,K)=B(L,K)+C(N)*TEMP
      GO TO 174
 702  J2=.FALSE.
      DO 271 K=1,NV
      TEMP=BE(1,K)
      RES=(FJ(K)-F1(K))/S
      BE(1,K)=RES
 271  B(1,K)=B(1,K)+(RES-TEMP)*W(1)
 174  CONTINUE
      IF(M.LT.NI) GO TO 175
      HSUM=ZERO
      VAL=TVAL**(-KE)
      DO 635 K=1,NV
      BDUBL=B(KE,K)
 635  HSUM=HSUM+BDUBL**2
      HSUM=VAL*DSQRT(HSUM)
      IF(NSF) GO TO 175
      IF(DABS((HSUM-SM)/HSUM).LT.0.01D0) GO TO 176
      SM=HSUM
 175  CONTINUE
!  THIS NEXT PART THE PROPER STARTING VALUE FOR T
 176  IF(NSF) GO TO 180
      TP=(SS/HSUM)**PW*DIR
      IF(NES) TP=XL
      IF(NES) GO TO 170
      IF(TP/T.GT.ONE) GO TO 170
      TP=0.8D0*TP
      NCOUNT=NCOUNT+1
      IF(NCOUNT.GT.10) RETURN
      GO TO 4000
 170  NSF=.TRUE.
!  FIND POSITION (AND VELOCITY FOR CLASS II AND IIS) AT THE END OF THE S
 180  DO 35 K=1,NV
      RES=B(KE,K)
      DO 34 L=1,KD
  34  RES=RES+B(L,K)
      X(K)=X(K)+V(K)*T+T2*(F1(K)*W1+RES)
      IF(NCL) GO TO 35
      RES=B(KE,K)*U(KE)
      DO 33 L=1,KD
  33  RES=RES+B(L,K)*U(L)
      V(K)=V(K)+T*(F1(K)+RES)
  35  CONTINUE
      TM=TM+T
      NS=NS+1
  74  IF(NPER) RETURN
!  CONTROL ON SIZE OF NEXT SEQUENCE AND RETURN WHEN TF IS REACHED
      GO TO 2222
      ENTRY CONT27
      NPER=.FALSE.
 2222 CALL FORCE(X,V,TM,F1)
      NF=NF+1
      IF(NES) GO TO 341
      TP=((SS/HSUM)**PW)*DIR
      IF(TP/T.GT.SR) TP=SR*T
 341  IF(NES) TP=XL
      IF (DIR*(TF-TM).GT.DIR*TP*1.001D0) GO TO 77
      TP=TF-TM
      NPER=.TRUE.
!  PREDICT B-VALUES FOR NEXT SEQUENCE,
  77  Q=TP/T
      DO 39 K=1,NV
      RES=ONE
      DO 39 J=1,KE
      IF(NS.GT.1) BT(J,K)=B(J,K)-BT(J,K)
      IF(J.EQ.KE) GO TO 740
      M=MC(J)
      JD=J+1
      DO 40 L=JD,KE
      B(J,K)=B(J,K)+XI(M)*B(L,K)
  40  M=M+1
 740  RES=RES*Q
      TEMP=RES*B(J,K)
      B(J,K)=TEMP+BT(J,K)
  39  BT(J,K)=TEMP
      NL=NI
      GO TO 722
      END
